//
//  Simulation.cpp
//  phy_four
//
//  Created by Snir Gazit on 3/12/12.
//  Copyright (c) 2012 __MyCompanyName__. All rights reserved.
//

#include "Simulation.h"

typedef boost::chrono::system_clock::time_point TimePoint;
typedef boost::chrono::duration<double, boost::ratio<3600>> Hour;
typedef boost::chrono::duration<double, boost::micro> MicroSec;

using namespace std;

double calcWR(double mu, double g, int m) {
	double c1, c2, c3, md = double(m), arg;
	double fac = 0.25 * (double(NF) + 2. * double(md));
	double fac2 = 0.25 * (2. + double(NF) + 2. * double(md));
	arg = 0.25 * pow(mu, 2) / g;
	c1 = 0.25 * pow(g, -fac2);
	c2 = mu * gsl_sf_gamma(fac2) * gsl_sf_hyperg_1F1(fac2, 3. / 2., arg);
	c3 = sqrt(g) * gsl_sf_gamma(fac) * gsl_sf_hyperg_1F1(fac, 1. / 2., arg);
	return c1 * (c2 + c3);
}

void calcWeight(double muIn, double gIn, int maxBond, ostream &outFile,
		BoltzWeightType &Rweight) {
	double mu = muIn, g = gIn;
	double n, c, w;
	c = calcWR(mu, g, 0);
	for (int i = 1; i < MAX_BOND; i++) {
		n = calcWR(mu, g, i);
		w = (n / c);
		c = n;
		Rweight[i - 1] = w;
	}
}

double calcWRGauss(double mu, int m) {
	double fac = double(NF) / 2 + m;
	return pow(mu, -fac) * gsl_sf_gamma(fac);
}

void calcWeightGauss(double muIn, int maxBond, ostream &outFile,
		BoltzWeightType &Rweight) {
	double mu = muIn ;
	double n, c, w;
	c = calcWRGauss(mu, 0);
	for (int i = 1; i < MAX_BOND; i++) {
		n = calcWRGauss(mu, i);
		w = (n / c);
		c = n;
		Rweight[i - 1] = w;
	}
}
// Worm Alg

void OnModel::init(QmcParameters &qmcParameters) {
	//Init Params
	qmcParameters_ = qmcParameters;
	string runNum = qmcParameters_.runNum.getVal();
	// OutFile
	string fileName = qmcParameters_.outPath.getVal() + "/out" + runNum
			+ ".txt";
	outFile_.open(fileName.c_str());
	outFile_ << "Starting N=3 Simulation, Run Number - " + runNum;
	// Init Lattice
	boost::array<int, DIM> sizes = qmcParameters_.Llist.getVal();
	lat_.init(sizes);
	lat_.isClosed_ = true;
	// Init Rng
	eng_.seed(qmcParameters_.seed.getVal());
	// Init Moves
	//	createMoveXY_.init(&eng_, &lat_, qmcParameters_, XYWeight);
	//	createMoveIsing_.init(&eng_, &lat_, qmcParameters_, IsingWeight);
	//  removeMoveXY_.init(&eng_, &lat_, qmcParameters_, XYWeight);
	//removeMoveIsing_.init(&eng_, &lat_, qmcParameters_, IsingWeight);
	createPairMove_.init(&eng_, &lat_, qmcParameters_, boltzWeight_);
	shiftMove_.init(&eng_, &lat_, qmcParameters_, boltzWeight_);
	chooseLoopType_.init(&eng_, UniIntVariableParams(0, NF - 1));
	setCoupeling(qmcParameters_.g.getVal());
	// Measurements
	measures_.init(&lat_, qmcParameters_, &boltzWeight_, &outFile_);
	// Read From File
	if (qmcParameters_.initFromFile.getVal())
		this->read();
	int totalNumberOfMeasure = (int) pow(
			(double) qmcParameters.binSize.getVal(),
			qmcParameters.numOfBins.getVal());
	//measures_.updateMinNumOfMeasure();
}
OnModel::~OnModel() {
	outFile_.close();
}

void OnModel::thermalize() {
	int numOfThermal = qmcParameters_.numOfThermal.getVal();
	acc CProb, XYProb;
	for (int i = 0; i < numOfThermal; i++) {
		nextMove();
		if (lat_.isClosed_)
			CProb(1.);
		else
			CProb(0.);
	}
	//CreateLoopMove_.setAvWormLength(2 * ba::mean(wormLength));
	double pClosed = ba::mean(CProb);
	outFile_ << "Final Closed Prob - " << pClosed << endl;
	int closedRate = int(double(pClosed) * double(lat_.getNumOfSites()));
	int openRate = int(
			double(qmcParameters_.rateFactor.getVal())
					* double(lat_.getNumOfSites()));
	//	closedRate = 1;
	measures_.updateClosedRate(closedRate);
	measures_.updateOpenRate(closedRate);
	outFile_ << "Closed Rate - " << closedRate << endl;
}

bool OnModel::run() {
	bool endTime, dumpFlag;
	namespace ba = boost::accumulators;
	ba::accumulator_set<double, ba::stats<ba::tag::mean> > wormLength;
	while (measures_.checkNotDone()) {
		nextMove();
		measures_.DoMeasurement(endTime, dumpFlag);
		if (dumpFlag) {
			this->dump(false);
		}
		if (endTime) {
			//this->dump(true);
			break;
		}
	}
	return !endTime;
}
void OnModel::nextMove() {
	static int wormLength = 0;
	if (lat_.isClosed_ == false) {
		if (shiftMove_.makeMove())
			wormLength++;
		if (lat_.isClosed_) {
			wormLength_[lat_.sectorType_](wormLength);
			wormLength = 0;
		}
	} else {
		int sectorType = chooseLoopType_.draw();
		createPairMove_.makeMove(sectorType);
	}
}
void OnModel::setCoupeling(double g) {
	if (qmcParameters_.simType.getVal() == "HEISENBERG")
		boltzWeight_.fill(1.);
	else if (qmcParameters_.simType.getVal() == "GAUSS")
		calcWeightGauss(qmcParameters_.mu.getVal(),
				qmcParameters_.maximalBonds.getVal(), outFile_, boltzWeight_);
	else
		calcWeight(qmcParameters_.mu.getVal(),g, qmcParameters_.maximalBonds.getVal(),
				outFile_, boltzWeight_);
	shiftMove_.setCoupeling(boltzWeight_);
	createPairMove_.setCoupeling(boltzWeight_);
}

void OnModel::simulatedAnnealing() {
	outFile_ << " Startting Simulated Anneling " << endl;
	double gHighT = qmcParameters_.gHighT.getVal();
	int numOfSMA = qmcParameters_.numOfGSA.getVal();
	double dg = (qmcParameters_.g.getVal() - gHighT) / double(numOfSMA - 1);
	double nLengthIsing, nLengthXY, nTotal;
	int numOfThermal = qmcParameters_.numOfThermal.getVal();
	vector<double> gSA;
	for (int i = 0; i < qmcParameters_.numOfGSA.getVal(); i++)
		gSA.push_back(gHighT + dg * double(i));
	for (auto it_g = gSA.begin(); it_g != gSA.end();) {
		double cg = (*it_g);
		setCoupeling(cg);
		outFile_ << "Current g -" << cg << endl;
		bool conv;
		int numIter = 0;
		for (auto it = wormLength_.begin(); it != wormLength_.end(); it++)
			*it = acc();
		int i = 0;
		while (i < numOfThermal) {
			nextMove();
				i++;
		}
		int tot = 0;
		for (int i = 0; i < NF; i++) {
			outFile_ << "N = " << i << " - " << ba::mean(wormLength_[i]) << " ";
			tot += ba::mean(wormLength_[i]);
		}
		numOfThermal = int(500 * double(lat_.getNumOfSites()) / double(tot));
		outFile_ << "Worm Length MA - " << tot << " numOfThermal "
				<< numOfThermal << " Iter " << numIter << endl;
		numIter++;
		//&& (numIter < 50));
//		if (conv)
		it_g++;
		//	else {
		//	if (it_g != gSA.begin())
		//	it_g--;
		//	}
	}
	int finished = 0;
	acc avTot;
	for (int i = 0; i < 80; i++) {
		do {
			nextMove();
			finished = count_if(wormLength_.begin(), wormLength_.end(),
					[](acc &a) {return (ba::count(a)>150);});
		} while (finished < NF);
		double tot =
				(1. / double(NF))
						* accumulate(wormLength_.begin(), wormLength_.end(), 0.,
								[](double a,acc&b)->double {return (a+ba::mean(b));});
		if (i > 60)
			avTot(tot);
		for (int j = 0; j < NF; j++) {
			outFile_ << "N = " << j << " - " << ba::mean(wormLength_[j]);
			tot += ba::mean(wormLength_[j]);
			wormLength_[j] = acc();
		}
		outFile_ << endl;
	}
	double tot = ba::mean(avTot);
	int closedRate = int(
			double(qmcParameters_.rateFactor.getVal())
					* double(lat_.getNumOfSites()) / tot);
	measures_.updateClosedRate(closedRate);
	measures_.updateOpenRate(lat_.getNumOfSites());
	outFile_ << endl << "Closed Rate - " << closedRate << endl;
	outFile_ << endl << "Finished Simulated Annealing " << closedRate << endl;
}

void OnModel::dump(bool final) {
	string outPath = qmcParameters_.outPath.getVal();
	measures_.dump(outPath, final);
	lat_.dump(outPath, final, &outFile_);
	eng_.dump(outPath);
	string fileName = outPath + "/sim_params.xml";
// Write to XML
	pugi::xml_document doc;
	pugi::xml_node rootNode = doc.append_child("SIM_PARAMS");
// Parameters
	pugi::xml_node node = rootNode.append_child("CLOSED_RATE");
	node.append_child(pugi::node_pcdata).set_value(
			(boost::lexical_cast < string > (measures_.getClosedRate())).c_str());
	doc.save_file(fileName.c_str());
//	doc.save(std::cout);
//	cout << endl;
}

void OnModel::read() {
	outFile_ << "Starting Reading From File" << endl;
	string outPath = qmcParameters_.outPath.getVal();
	measures_.read(outPath);
	outFile_ << "Measures Loaded" << endl;
	lat_.read(outPath);
	int totalBonds;
	try {
		bool testConf = lat_.testConfiguration(totalBonds);
		outFile_ << "Total Bonds Calculated - " << totalBonds
				<< " total Bonds previous sim - " << lat_.totalBonds_ << endl;
		if (totalBonds != lat_.totalBonds_)
			throw "Total Bonds after Read don't match";
		if (!testConf)
			throw "Bad Load of Configuration";
	} catch (string &e) {
		outFile_ << e;
		exit(1);
	}
	outFile_ << "Lattice Loaded" << endl;
	eng_.read(outPath);
	string fileName = outPath + "/sim_params.xml";
// Read From XML
	pugi::xml_document doc;
	try {
		pugi::xml_parse_result result = doc.load_file(fileName.c_str());
		if (!result)
			throw result;
	} catch (pugi::xml_parse_result &result) {
		std::cout << "Error description: " << result.description() << "\n";
		exit(1);
	}
	pugi::xml_node rootNode = doc.child("SIM_PARAMS");
	measures_.updateClosedRate(readXmlData<int>(rootNode, "CLOSED_RATE"));

}
